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ABSTRACT 

This  paper  investigates  the  spectrum  of  the  iteration  operator  of  some  finite  element  pre¬ 
conditioned  Fourier  collocation  schemes.  The  first  part  of  the  paper  analyses  one-dimensional 
elliptic  and  hyperbolic  model  problems  and  the  advection-diffusion  equation.  Analytical  ex¬ 
pressions  of  the  eigenvalues  are  obtained  with  use  of  symbolic  computation.  The  second 
part  of  the  paper  considers  the  set  of  one-dimensional  differential  equations  resulting  from 
Fourier  analysis  (in  the  tranverse  direction)  of  the  2-D  Stokes  problem.  All  results  agree  with 
previous  conclusions  on  the  numerical  efficiency  of  finite  element  preconditioning  schemes. 
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1.  INTRODUCTION 


In  the  recent  past,  Chebyshev  collocation  schemes  have  been  applied  extensively  to  the 
numerical  integration  of  the  Navier-Stokes  equations  [1,  3,  4].  For  scalar  elliptic  problems,  it 
is  well  known  that  the  condition  number  of  the  matrix  system  of  discrete  algebraic  equations 
increases  rapidly  with  N,  the  number  of  degrees  of  freedom  of  the  problem  at  hand.  Therefore, 
the  preconditioning  technique  seems  to  be  the  only  adequate  tool  in  order  to  overcome 
this  numerical  burden.  The  present  authors  [5,  6]  demonstrated  that  finite  elements  (FE) 
constitute  powerful  preconditioners  for  general  second-order  elliptic  equations.  In  [3],  several 
fluid  flow  elements  in  velocity-pressure  formulation  were  investigated.  From  the  analysis  of 
the  eigenspectrum  of  the  iteration  operator,  it  was  shown  that  the  Q2-Q1  element  is  the 
best  choice  for  the  steady  Stokes  problem.  As  all  the  previous  analyses  on  finite  element 
preconditioning  were  carried  out  numerically,  the  present  note  aims  at  analytical  results 
through  use  of  symbolic  manipulation  languages  (cfr.  [10]). 

For  the  sake  of  simplicity,  we  will  restrict  ourselves  to  the  study  of  a  finite  element 
preconditioned  Fourier  collocation  scheme.  In  this  case,  the  mesh  size  is  uniform  and  the 
algebra  is  considerably  reduced.  In  Section  2,  a  one-dimensional  model  is  considered.  The 
collocation  process  is  preconditioned  by  Lagrangian  linear,  quadratic,  cubic  and  Hermite 
cubic  elements,  respectively.  The  Richardson  iteration  method  is  set  up  with  these  FE 
preconditioners  as  approximate  operators  and  algebraic  solvers.  Using  the  spatial  structure 
of  the  eigenvectors  of  the  Fourier  solutions,  one  may  perform  a  full  analysis  of  the  eigenvalues 
of  the  iteration  operator.  This  theoretical  investigation  corroborates  the  previous  numerical 
analyses  [6].  In  Section  3,  a  one-dimensional  hyperbolic  model  is  investigated  using  linear  and 
quadratic  FE  preconditioning.  The  upwinding  technique  is  also  examined.  A  further  model 
consists  in  an  advection-diffusion  equation.  In  Section  4,  the  Stokes  problem  is  reduced  to 
a  1-D  incompressible  flow  model  amenable  to  Fourier  discretization.  The  Q2-Q1  and  Ql-PO 
elements  are  candidates  as  preconditioners.  A  similar  Fourier  analysis  is  done.  The  results 
corroborate  numerical  experiments  carried  out  in  the  framework  of  Chebyshev  collocation 
[3], 


2.  ELLIPTIC  MODEL 

l>t  us  first  consider  the  simple  elliptic  problem: 

-u«  =  /,  0  <  i  <  2tt,  (2.1) 

with  periodic  boundary  conditions.  The  subscript  indicates  partial  derivative.  The  Fourier 
approximation  of  the  dependent  variable  u  is 

N/l-\ 

uN=  £  upeiFI'»  0  <J<N,  (2.2) 

p=-N/2 

where  Up  are  the  discrete  Fourier  coefficients  and  Xj  the  collocation  points  defined  by 

xi  =  i  G  [0,  iV[.  (2.3) 
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The  linear  system  corresponding  to  (2.1)  may  be  found  in  [1]  and  will  be  denoted  by  Lc. 
The  eigenfunctions  of  (2.1)  are 

6(p)  =  e*%  0  <j<N,  (2.4) 

with  the  corresponding  eigenvalues 

X(p)=p\  pG[-y,y-l].  (2.5) 

The  collocation  problem  will  be  preconditioned  by  finite  elements.  Introducing  the  approxi¬ 
mate  FE  operator  L,  the  preconditioned  Richardson  iteration  is  written  as: 

uk+1  =  uk-  akL~\Lcuk  -  /),  (2.6) 


where  k  is  an  iteration  index,  ak  a  relaxation  factor  and  u,  /  the  vectors  corresponding  to  the 
unknowns  and  source  terms  at  the  collocation  points.  The  convergence  of  (2.6)  is  governed 
by  the  spectral  radius  p(A)  of  the  iteration  operator  defined  by  A  —  I  —  akL~1Lc.  The 
optimal  value  of  the  relaxation  factor  is: 

a°pt  =  7  7~ T  >  (2-7) 

Amin  i  ^moi 


where  Amm  and  Xmax  are  the  minimum  and  maximum  eigenvalues  of  L  1LC.  An  approximate 
estimate  of  the  number  of  iterations  n  needed  to  reduce  the  error  norm  by  a  factor  (  is  given 

by 

n  =  -log(/R00{A),  (2.8) 


where  Roo(-A)  =  —  logp(A)  is  the  asymptotic  rate  of  convergence  of  the  iteration  matrix. 
The  spectral  radius  p(A)  which  is  involved  in  the  error  reduction  process  with  the  use  of  a^t 
(Eq.  (2.7))  is  given  by 

p(A)  =  (2.9) 

Amax  i  Amin 

In  order  to  investigate  this  quantity  for  various  preconditioners,  we  have  to  define  the 
finite  element  problem  more  precisely. 

Lagrangian  linear  elements,  Hermite  cubic  elements  (i.e.,  Ql,  P3  in  Ciarlet’s  notations  [2]) 
as  well  as  higher-order  Lagrangian  interpolants  have  their  vertices  at  the  Fourier  collocation 
grid  (2.3).  However,  for  Lagrangian  quadratics  (Q 2),  mid-points  are  added  at 


27T..  1 

xit\  ~  A  ^  +  2  ’ 


j  6  [0,  JV[, 


(2.10) 


while  for  Lagrangian  cubics  (Q3),  we  have  additional  grid  nodes  located  at 


j  e  [0, JV[,  (2.11) 

with  obvious  definitions.  For  the  Lagrangian  case,  the  FE  unknowns  are  the  nodal  values, 
while  for  the  Hermite  case,  the  unknowns  are  u;  and  u',  the  prime  denoting  the  first-order 
derivative  of  u.  In  [5],  the  iteration  operator  is  written  as: 


A  =  1  -  S^MhlhLe, 


(2.12) 


2 


where  Sh  is  the  stiffness  matrix,  Mh  the  mass  matrix,  I h  an  interpolation  matrix  evaluating 
the  Fourier  interpolant  of  the  collocation  operator  at  the  FE  nodes.  The  mesh  size  of  the 
FE  grid  is  defined  by 

h  =  2  n/N. 


For  Ql  elements,  Ih  reduces  to  the  unit  matrix;  for  higher-order  interpolants  however,  the 
structure  of  this  matrix  is  more  complicated.  In  order  to  avoid  writing  the  details  of  Ih, 
we  will  systematically  assume  the  use  of  static  condensation.  Consequently,  the  iteration 


operator  may  be  written: 

A  —  I  -  Sh'MhL,, 


(2.13) 


where  Sh  and  Mh  refer  to  stiffness  and  mass  matrices  after  static  condensation. 


2.1.  Linear  Lagrangian  Elements 


For  an  interior  node,  the  expressions  for  the  stiffness  and  mass  matrices  are  well-known: 

ShU:  =  +  2 Uj  -  ui+1),  (2.14) 

MJ,  =  £(/>->+ 4/; +  W  (2.15) 

Fourier  analysis  of  (2.14),  (2.15)  with  the  eigenfunction  (2.4)  leads  to  the  expression  of  the 
eigensnectrum  of  S^1  MhLc,  denoted  by  cr(p), 


<r(p)  = 


(p/i/2)2  (2  +  cosph) 
sinJph/2  3 


N  ^  <rN  i 

- <  p  < - 1. 

2  2 


(2.16) 


Typically,  the  second  factor  in  the  right-hand  side  of  (2.16)  comes  from  the  contribution 
of  the  mass  matrix.  In  the  case  of  finite  difference  (FD)  preconditioning,  this  factor  is 
one.  For  p  =  0,  o(p)  =  1,  while  for  p  =  —N/2,a(p)  =  7r2/12.  This  last  value  should  be 
compared  to  the  FD  equivalent  which  is  cr(p)  =  7t3/4  [7).  The  eigenvalue  spectrum  of  the 
FE  preconditioning  is  reduced  because  of  the  beneficial  presence  of  the  mass  matrix.  Fig.  1 
shows  the  behavior  of  cr(p)  with  respect  to  p  for  h  —  27r / 100.  The  function  has  a  minimum 
value  equal  to  0.693.  Therefore,  the  optimum  value  for  a  is 


ctopt  s  1.18,  (2.17) 

and  over-relaxation  is  possible  for  FE  preconditioning  unlike  the  FD  preconditioning  where 
under-relaxation  is  required  to  converge.  In  practice,  the  Ql  preconditioning  with  a  spectral 
radius  of  0.18  converges  twice  as  fast  as  the  FD  preconditioner  whose  spectral  radius  is  of 
the  order  of  0.42. 
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2.2.  Quadratic  Lagrangian  Elements 


The  equations  related  to  nodes  j  and  j  ±  ~  may  be  cast  in  the  following  matrix  form: 

/  \ 

l,  (  1  8  10  0  \  |  /,_*  | 

(2.18) 
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The  use  of  static  condensation  eliminates  the  contribution  of  u^_i  and  uJ+i  and  Eq.  (2.18) 
reduces  to  only  one  relationship  for  node  j  on  the  collocation  gri^ 


+  2uj  ~  ui+i)  =  i  +  fi  +  i)- 


(2.19) 


Let  us  notice  that  for  Sh  on  the  left  hand-side  of  (2.19),  we  recover  the  stiffness  matrix  Sh 
associated  to  Q1  elements  whereas  in  the  right-hand  side,  Mh  corresponds  to  a  different 
quadrature  rule.  Carrying  out  the  Fourier  analysis  of  (2.19),  one  obtains 


CT(?,)  =  +  2coa(pft/2>)’ 


2  ~  y  ~  2 


(2.20) 


For  the  particular  vaiues  p  =  0  and  p  =  -N/2,a(p)  is  equal  to  1  and  7r2/ 12, respectively.  As 
o(p)  is  a  monotically  decreasing  function  with  respect  to  p  (Fig.  1),  the  optimum  value  of  a 
is 

aopt  —  2/(1  +  tt2/12)  =  1.0974,  (2.21) 

and  the  corresponding  spectral  radius  p(A )  is  equal  to  0.0974. 


2.3.  Cubic  Lagrangian  Elements 

For  the  sake  of  compactness,  we  give  the  local  stiffness  and  mass  matrices  over  the  uniform 
mesh: 


Sh  =  - 


1 

r  37/20 

-189/80 

27/40 

-13/80  \ 

I 

-189/80 

27/5 

-297/80 

27/40 

i 

27/740 

-297/80 

27/5 

-189/80 

s,  -13/80 

27/40 

-189/80 

37/20 

) 

/  16/105 

33/260 

-3/70 

19/840  \ 

h 

33/280 

27/35 

-27/280 

— 3/70 

2 

-3/770 

-27/280 

27/35 

33/280 

\  19/840 

-3/70 

33/280 

16/105/ 

(2.22) 


(2.23) 


Assembling  the  matrices  of  (2.22),  (2.23)  over  two  adjacent  elements  and  eliminating  the 
four  unknowns  attached  to  the  interior  nodes  Uj±\/3,Uj±2/3,  we  are  left  with  the  relation: 

+  2u,  -  uitl)  =  A( III  + 1/,_, + 3/y_,  +  f /, + s/it,  +  |/i+l  +  f-f).  (2.24) 
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Here  again,  as  in  the  previous  case,  we  recover  in  the  left-hand  side  of  (2.24)  the  stiffness 
matrix  of  Ql  elements.  Fourier  analysis  of  (2.24)  leads  to  the  eigenvalue  spectrum  a(p)  . 


a(p)  = 


(ph/2)3  1.4  3  ph  ,  vh  ph  2. 

SW)Io(73CO!  T  +  3cos  t+5“sT  +  3)' 


N  N 

-y  <P<  y-1.  (2.25) 


The  particular  values  of  a(p)  corresponding  to  p  =  0  and  p  =  —N/2  are  cr(p)  =  1  and 
497r2/480  ~  1.007522.  The  optimal  value,  of  the  relaxation  factor  is  given  by  2/(1  + 
497r2/480)  ~  0.9963  and  the  corresponding  spectral  radius  p(A )  is  equal  to  0.00375. 

Looking  back  at  the  results  in  the  previous  subsections,  one  observes  that  the  spectral 
radius  p( A)  diminishes  with  increasing  polynomial  degrees.  This  does  not  mean  however 
that  one  should  use  higher  order  elements  in  the  preconditioning  because  they  involve  more 
computational  work  as  the  bandwidth  of  the  algebraic  system  increases. 


2.4.  Cubic  Hermite  Elements 


At  node  j,  the  discrete  equations  are 


h  9  26  9 

+  7(Io/j-1  +  T/2  +  io/j+l)’ 

^j(«j-i  -  ui+i)  ~  ~  H  +  u',+i )  =  ~  /i+i) 


30' 


420 

+  /;,:)• 


Fourier  analyzing  (2.26),  one  gets  the  spectrum: 

(ph/2)3 


(2.26) 


(2.27) 


j(p)  = 


6  sin2(ph/2)  —  (ph/2)  sin(ph/2)  cos(ph/2)  7 


1  13 

-(26  +  9  cos  ph  +  —  ph  sinp/i), 


\ph\  €  [0,4  (2.28) 

For  p  =  0  and  —N/2,  a(p)  =  1  an  177r2/168  ~  0.9987,  respectively.  Fig.  1  displays  the 
behavior  of  cr(p),  which  is  first  slightly  decreasing  with  respect  to  p  achieving  its  minimum 
value  0.97722  and  then  increasing  to  1.  The  optimal  value  of  a  is  1.01152  and  the  cor¬ 
responding  spectral  radius  0.0115,  an  intermediate  value  in  between  those  of  Q2  and  Q3 
preconditioning. 


3.  HYPERBOLIC  PROBLEMS 

We  now  turn  our  attention  to  the  first-order  differential  equation 

u*  =  /. 


(3.1) 
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in  the  periodic  case.  The  eigenfunctions  of  (3.1)  are  again  (2.4)  with  the  eigenvalues 

a(p)  =  *p>  p  e  l-y -  !]• 


3.1.  Linear  Lagrangian  Elements 

Using  the  standard  Galerkin  approach,  a  centered  scheme  is  produced  and  yields  the 
discrete  equation 

+  4  f,  +  /,«).  (3.2) 


By  Fourier  analysis,  we  obtain: 


ph  2  +  cos  ph 
sinp/i  3  ’ 


\ph\  e  [0,7 r]. 


For  p  =  0,  cr(p)  =  1,  while  a(—N/2)  is  obviously  unbounded  as  in  the  FD  case.  The  presence 
of  the  mass  matrix  does  not  help  to  circumvent  the  difficulty. 

One  may  proceed,  however,  using  an  upwinding  technique.  This  has  been  a  key  step 
to  treat  hyperbolic  problems.  In  finite  elements,  the  method  uses  separate  test  and  trial 
function  spaces,  i.e.,  the  Petrov-Galerkin  method.  There  are  several  ways  to  implement 
upwinding.  Let  us  introduce  the  weight  functions  U7*(r),(i  =  1,2)  defined  on  the  reference 
interval  [-1,+1]  by  Heinrich  and  Zienkiewicz  [8]: 


u>‘(r)  =  <p‘(r)  +  (  — l)‘eF(r), 


-1  <  r  <  1, 


where  <p‘(r)  are  the  linear  Lagrangian  trial  functions,  F(r )  an  auxiliary  quadratic  element 
vanishing  at  both  end  points 

F(r)  =  l(l-r1).  (3,5) 

and  e  the  upwinding  parameter  to  be  given  independently.  Using  <p‘  and  w'  as  trial  and  test 
functions  respectively,  one  gets 


+  iUi + 


-  2  “f*  3C  _  Lt  -  4a  -  WV, 

'Uj+i  =  ~ w f,~i  +  3/j  +  ~Trfj+i- 


2  —  3e 
12  J 


This  equation  reduces  to  (3.2)  when  e  =  0  (no  upwinding).  With  a  value  of  e  as  yet  undefined, 
the  Fourier  analysis  of  (3.6)  leads  to  complex  eigenvalues  given  by: 


ff(P)  =  C 


6  e(  1  —  cos  ph)  +  i  sin  ph 


(3eph  sinph  +  2iph(2  +  cosph)),  |pA|  E  [0,7r].  (3.7) 


For  p  =  0,  a(p)  =  1  while  for  p  =  — iV/2,a(p)  =  i/2  independently  of  the  value  of  e. 
Upwinding  forces  the  spectrum  of  S^1  MhLc  almost  entirely  inside  the  unit  circle,  as  shown 
on  Fig.  2  where  the  eigenvalues  (3.7)  have  been  plotted  for  e  =  1.  The  spectral  radius 
of  the  matrix  A  in  this  case  is  equal  to  \/b/2  and  under-relaxation  is  required  in  order  to 
ensure  convergence  of  the  preconditioning  iterations.  The  eigenvalues  (3.7)  being  complex, 
the  evaluations  of  a  opt  and  the  spectral  radius  p(A )  are  no  longer  given  by  (2.7)  and  (2.9). 
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3.2.  Quadratic  Lagrangian  Elements 


Another  way  to  solve  (3.1)  consists  in  using  a  FE  preconditioner  based  on  quadratic 
Lagrangian  elements.  Applying  the  Galerkin  approach  and  assembling  the  contributions  of 
two  adjacent  elements  at  node  j ,  one  obtains  a  set  of  three  equations  related  to  nodes  j  and 
j  ±  5  similar  Lo  (2.18),  which  are  cast  in  matrix  form: 


(  u,~l  \ 

(  fj-i  \ 

0  0  \ 

<1  1  1 

Ui  l 
•'“I 

h 

/  2  16  2  0  0  \ 

3  ~  6 

Uj 

~  30 

—  1  2  8  2  —1 

fi 

0  I  ) 

V  0  0  2  16  2  / 

'  f]+i  ' 

(3.8) 


If  static  condensation  is  carried  out  through  Gaussian  elimination,  the  contribution  of  the 
exterior  nodes  Uj-i,itj+i  disappears  and  one  is  left  with  a  staggered  scheme: 


~  ui+j)  -  +  124-j  +  18/>  +  12/j+i 


-  fi+ 1)> 


for  Fourier  analysis.  Its  spectrum  is  given  by: 


a(p)  = 


f  9  +  12  cos  ^  —  cos  ph 
sinf  20 


\Ph\  e  [o,7r]. 


(3.9) 


(3.10) 


It  is  monotonically  decreasing  and  bounded  by  a(0)  =  1  and  a(  —  N/ 2)  =  7r/4  as  shown  on 
Fig.  3.  One  should  notice  that  the  first  term  in  the  right-hand  side  of  (3.10)  is  identical 
to  the  spectrum  obtained  in  the  FD  case  where  the  function  is  computed  on  the  main  grid 
while  the  derivative  is  evaluated  by  first-order  differences  on  a  staggered  grid.  The  second 
term  whose  maximum  value  is  equal  to  1  is  induced  by  the  presence  of  the  mass  matrix  and 
reduces  to  unity  in  the  FD  case.  The  optimal  value  of  a  is  equal  to 


<*opt  =  2/(1  +tt/4)  =  1.1202, 

and  the  corresponding  spectral  radius  p(A)  is  equal  to  0.1202.  This  staggered  scheme  gener¬ 
ated  by  Q2  elements  is  the  key  of  success  for  FE  preconditioning  of  Navier-Stokes  problems. 
This  excellent  behavior  explains  the  reason  why  in  Demaret-Deville  [5],  the  relaxation  pa¬ 
rameter  was  almost  independent  of  the  Reynolds  number. 


3.3.  Advection-Diffusion  Model 

The  last  scalar  model  analyzed  in  this  paper  is  the  one- dimensional  advection-diffusion 
problem.  The  differential  equation  writes 

-kuxx  +  cux  =  f(x),  (3-11) 

where  k  is  the  diffusion  coefficient  and  c  the  constant  advection  velocity.  Particular  interest 
bears  on  advection  dominated  problems  which  impose  severe  conditions  on  the  element  mesh 
size  (cfr.  Thomasset  [9]).  With  the  eigenvectors  (2.4),  the  eigenvalues  of  (3.11)  are 


where  7  is  the  cell  Reynolds  number  defined  by  7  =  chjt c. 

Using  the  linear  FE  basis  and  upwinding  introduced  in  the  previous  section,  one  gets  the 
discrete  equations 

~(1  +  )«>-i  +  (2  +  -ye)u}  -  (1  -  7“"y~)ui+i  =  ^[(2  +  3  <:)/,_  1  +  8/,  +  (2  -  3e)/i+1]. 

(3.12) 

where  e  is  the  upwinding  parameter  of  Eq.  (3.4).  With  no  upwinding  (i.e.,  e  =  0),  stability 
requirements  restrict  7  to  values  <  2.  The  Fourier  analysis  of  (3.12)  is  straightforward.  The 
eigenvalues  of  (3.12)  are  complex  and  given  by: 


.  ^^-(2  +  cospA)  +  £7^  sinph  +  i [7 ^ ( 2  +  cos ph)  -  sinp/i] 

P  2(2  +  erf)  sin2  ^  +  17  sin  ph  ’ 

In  absence  of  upwinding,  Eq.  (3.13)  reduces  to  an  analytical  expression 
imaginary  parts  may  be  written  in  compact  form: 


\ph\  e  [o,7r], 

(3.13) 
whose  real  and 


Re(cr(p))  =  ph 


Aph  sin2  ^  +  72  sin  ph  2  +  cos  ph 
16  sin4  ^  +  72  sinp/i  3 


Im(o{p))  —  'yph 


4  sin2  —  ph  sin  ph  2  -f  cos  ph 
16  sin4  ^  +  72  sin  ph  3 


(3.14) 


The  factor  (2  +  cos  ph)/ 3  in  the  right-hand  side  of  (3.14)  is  another  example  of  the  contri¬ 
bution  of  the  mass  matrix  in  FE  preconditioning.  Like  in  Eq  (2.16),  this  factor  reduces 
to  unity  in  the  expression  of  the  eigenvalues  corresponding  to  FD  preconditioning.  One 
can  draw  similar  conclusions  to  the  diffusion  problem,  except  for  the  complex  nature  of  the 
eigenvalues.  Introducing  ph  —  0  and  ph  —  n  into  the  eigenvalues  of  the  FD  case  gives  the 
bounds  of  the  spectrum: 


1  <  Re(aFD)  <  0  <  Im(aFD)  <  (3.15) 

In  the  FE  case,  the  upper  bounds  are  reduced  by  a  factor  3  because  of  the  presence  of  the 
mass  matrix.  Figure  4  displays  the  result  (o.l4)  for  both  FD  and  FE  preconditionings  and 
for  two  different  values  of  7  (i.e.,  0.2  and  2).  Even  in  the  limit  case  7  =  2,  the  spectrum  of 
A  for  FE  preconditioning  lies  inside  the  unit  circle.  Reducing  the  value  of  the  cell  Reynolds 
number  brings  the  eigenvalues  closer  to  the  real  axis. 

Figure  5  exhibits  (3.13)  with  7  =  2,  with  and  without  upwinding.  In  the  upwinding 
case,  €  was  chosen  equal  to  1.  The  eigenspectrum  is  rotated  counterclockwise  and  slightly 
stretched  inducing  a  somewhat  larger  spectral  radius. 


4.  STOKES  EQUATIONS 

Let  us  write  the  Stokes  equations  in  stress  formulation: 

divg.  + pf_  =  0,  (4.1) 
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divv  =  0. 


The  symbol  a  denotes  the  stress  tensor,  p  is  the  density,  /  the  body  force  term  and  v  is 
the  velocity  field.  Eq.  (4.1)  is  the  momentum  equation  and  Eq.  (4.2)  enforces  the  continuity 
constraint.  The  2-D  Stokes  problem  may  be  reduced  to  a  1-D  problem  if  the  solution  of 
(4.1),  (4.2)  is  sought  as  a  Fourier  mode: 


v  —  y(x)eikv,p  =  p(x)e'ky.  (4.3) 

Introducing  (4.3)  in  (4.1),  (4.2),  we  get: 

dp  du 2  dv 

d  dv 

p—(iku  +  -^)  -  ikp  -  2 pk2v  +  pfv  —  0,  (4.5) 

du 

— — h  ikv  =  0,  0  <  x  <  1-k.  (4.6) 

Ox 

The  velocity  and  pressure  fields  are  assumed  to  be  27r-periodic.  This  1-D  problem  is  dis¬ 
cretized  in  the  x  direction  using  Fourier  series  of  type  (2  2)  for  each  variable.  The  discrete 
collocation  equations  are  preconditioned  by  finite  elements  such  as  the  Q 2  —  Q 1  and  Ql  -  P0 
elements.  The  FE  equations  come  from  Galerkin  projection.  Introducing  ipi  and  <p;  the  trial 
functions  for  the  FE  approximations  of  the  velocity  and  pressure  fields,  respectively,  such 
that 

M  N 

n(*)  =  HlhV'i,  p(z)  =  £pi<Pi,  (4-7) 

1=1  1=1 

the  discrete  FE  equations  are  obtained  by  use  of  the  divergence  theorem  as  tool  for  the 
integration  by  parts  with  the  notation  fx  —  f ,  /„  =  g: 

+  -  ikpJ2Cj‘v‘  ~  DJ‘P‘  =  0  <  J  <  Af,  (4.8) 

i  ill 

ikp^CjfUi  +  ^[2Mfc25j(  +  pAji]vt  +  ik^  Ejipi  =  B]tgi,  0  <  j  <  M,  (4.9) 


-  H  DJiu‘  ~lkYl  EJiv<  =  °>  0<j<N. 

i  i 

In  (4.8-4.10),  the  various  matrices  are  defined  by  the  relationships: 

AP  =  J  ~^7~^dx>  =  J  =  J 

dj,= j -fcdx'  Eii = /***. 


(4.10) 


(4.11) 
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4.1.  Q2-Q1  Elements 


For  this  element,  M  =  2 N.  Carrying  through  the  algebra  involved  by  the  quadratures 
(4.11)  and  assembling  by  direct  stiffness  the  contributions  of  the  two  elements  connected  to 
node  j,  we  obtain: 


2u.  . 

~  i/2  +  14 Uj  —  8uj+ 1/2  4-  Uj+i) 

+  ~30~~(  — ' 0,-1  +  2Uj_1/2  +  +  U7  +  l)  -  -  4Vj_i/2  +  4v,>i/2  -  UJ  +  i) 


-g  (Pi-1  -  Pj+i)  -  ^q(-/j-i  +  2/j_i/2  +  8/j  +  2/j+i/2  -  /j+i), 

4u  uk2h.  .  2ik^i  . 

— (-4u;  +  8u;+1/2  -  4u,-+i)  +  -^-(2^  +  16ui+1/j  +  2uJ+i) y- (-v>  +  vi+V 

2  /i 

~  g(Pj  -  Pj+i)  =  20 (2/;  +  16/>+i/2  +  2fj+i), 


(4.12) 

(4.13) 


ikfi 


(~uj- 1  +  4tij_ 


1/2 


fik2h 

4uj  +  j/2  +  Uj  +  i) - jg-(- 


-U;_1  +  2u_,_i/2  +  +  2Vj+i/2  —  Vj+i 


2/1/  n  \ 

+  y(Vi  -  8uj-i/2  +  14u;  ~  8u/+i/2  +  v,-+i)  -  -yP7 

=  2g(-J>-i  +  2 +  8g7  +  2<7J+1/3  -  p>+i), 

2  ik^x  .  8/x,  ,  2  pk2h  v 

— 2~ (u>  “  U7  +  0  +  +  2uJ  +  1/2  ~  U7  +  l) jy- (U7  +  8v>+1/2  +  Vj+l) 

-~(Pj  +  Pjt-l)  =  -(/;  +  8/;+i/2  +  /j+l)r 

1  2ifc  . 

-(-Uj_i  -  4uj_i/2  4-  4uj+1/2  +  Uj+i)  4-  —  (v>-i/a  4-  Vj  +  Uj+i/a)  = 


(4.14) 

(4.15) 


(4.16) 


Eqs.  (4.12),  (4.13)  and  (4.16)  correspond  to  momentum  and  incompressibility  relations, 
while  (4.13)  and  (4,15)  are  the  momentum  equations  associated  to  mid-node  xy+i/2-  Similar 
expressions  hold  for  mid-node  x7_i/2  with  appropriate  shifts  for  the  indices. 

Now,  static  condensation  represents  a  formidable  task  and  is  greatly  helped  by  the  sym¬ 
bolic  manipulation  program.  Elimination  of  u7±i/a,  vj±i/2  leads  to  a  matrix  system  of  order 
3. 

The  full  Fourier  rolution  gives  the  collocation  matrix  Le\ 


Lc 


1  —  2fil 2  —  k2fi  —pkl  — il ^ 

—fikl  ~nl2  —  2fik 2  —ik 
\  -il  -ik  0  ) 


(4.17) 
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where  l  is  the  wavenumber  in  the  i  direction. 

The  analytical  computation  of  S^lMhLc  is  performed  as  far  as  the  symbolic  program 
cam  handle  tractable  expressions.  Then  numerical  evaluation  of  the  eigenspectrum  is  done. 
Because  of  the  divergence-free  constraint,  a  zero  eigenvalue  is  systematically  obtained.  In 
Figures  6  and  7,  the  eigenspectrum  of  S^1  MhLc  are  plotted  for  two  cases  k  =  1  and  k  = 
10,  respectively.  In  these  two  figures,  the  lower  curve  shows  the  same  behavior  as  the 
eigenspectrum  of  the  elliptic  problem  solved  by  Q2  elements.  For  l  =  -N/2,  a(l)  is  equal  to 
vr3/12.  For  the  other  curve,  cr(0)  =  1  and  a(—N/2)  is  close  to  2.07.  Therefore,  the  optimal 
a  value  is  given  by 

Q-opt  —  2/(2. 07  +  7r2 /1 2)  ~  0.69, 

a  value  close  to  2/3  obtained  by  Demaret-Deville  [3]  for  a  2-D  Chebyshev  collocation  dis¬ 
cretization  of  the  Stokes  problem  preconditioned  by  Q2-Q1  elements. 


4.2.  Q1-P0  Element 


The  quadratures  (4.11)  provide  less  complicated  discrete  equations  in  this  case: 


2jj.  fihk2  ika. 

y  (-«>-i  +  2 u,  -  Uj+i) - —  (ttj_i  +  4m,  +  uJ+1)  +  —  (vJ+1  -  Vj. x) 

-Pj  +  x  =  ^(/i- I  +  4/;  +  f j+ 1 ) ) 


(4.18) 


ikp  g.k2h.  n,  . 

—  (“i-i  -  u:+i) 3 — (v,-_i  +  4uj  +  vi+l)  +  +  2vj  ~  u>+i) 

ikh.  ,  h 

t  Pi+i)  -  g(3j-i  +  4 g,  +  5j+i ), 


(4.19) 


ikh 

~uj-i  +  u.7+i  +  y- (^-i  4"  2vJ  4“  vj+1)  ~  O'  (4.20) 

Obviously,  this  element  generates  second-order  differences  for  partial  derivatives.  When  the 
mass  matrix  is  involved,  the  standard  weighted  mean  between  three  adjacent  nodes  appears 
in  the  expressions.  No  static  condensation  is  needed  in  this  case.  Fourier  analyzing  Eqs. 
(4. 18)-(4.20),  the  stiffness  and  mass  matrices  are  now: 


Sh  = 


in  sin2(y)  +  amp  +  cos  hi) 
j. ik  sin  hi 
\  2i  sin  hi 


[ik  sin  hi 

+  cos  hi) 
ikh(  1  4-  cos  hi) 


2isin(y)  \ 
ikh cos(j)  I  i 
0  / 


Mh  =  diag(—( 2  +  cos  hi),  ^(2  +  cos  hi),  0).  (4-21) 

u  O 

In  Figures  8  and  9,  the  eigenspectrum  of  S^lMhLc  are  displayed  for  k  =  1  and  10,  respec¬ 
tively.  In  these  two  figures,  the  top  curve  is  that  of  the  elliptic  model  preconditioned  by  Ql 
element.  The  bottom  curve  starts  from  1  for  /  =  0,  decreases  till  a  minimum  value  close  to 
0.49  and  then  increases  to  reach  cr(  — IV/ 2)  =  0.5.  The  optimum  value  is 


aopt  =  2/(1  +  0.5)  =  4/3. 
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5.  CONCLUSIONS 


In  this  paper,  we  have  Fourier  analyzed  the  eigenspectrum  of  the  iteration  operator  for 
finite  element  preconditioning  of  Fourier  collocation  applied  to  one-dimensional  problems. 
For  elliptic  models,  this  theoretical  analysis  confirms  previous  numerical  findings,  especially 
the  beneficial  presence  of  the  mass  matrix  which  reduces  the  bounds  of  the  eigenspectrum. 
For  first-order  problems,  linear  elements  without  and  with  upwinding  are  considered.  With 
quadratic  elements,  a  staggered  scheme  is  produced.  Its  eigenspectrum  is  bounded  and  ranges 
between  1  and  tt/4.  Finally,  a  Stokes  problem  is  reduced  to  a  one-dimensional  approach. 
Two  types  of  elements  are  examined.  The  Q2-Q1  element  leads  to  an  optimum  value  of 
the  relaxation  parameter  close  to  the  value  obtained  by  numerical  analysis  of  preconditioned 
Chebyshev  collocation.  For  Ql-PO  element,  the  method  can  be  over-relaxed. 
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Figure  1.  Eigenspectrum  of  $hlMhLc  for  an  elliptic  model.  Preconditioners:  Ql:*;  Q2: 
Q3:A;  Cubic  Hermite:  Q. 
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Figure  2. 
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Figure  4.  Eigenspectrum  of  an  advection-diffusion  problem.  The  curves  outside  the  unit 
circle  are  concerned  with  FD  preconditioning  while  the  inside  curves  are  related  to  FE 
preconditioning.  Both  top  curves  are  obtained  for  the  cell  Reynolds  number  7  =  2.  The 
bottom  curves  are  gotten  for  7  =  0.2. 
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Figure  5.  Eigenepectrum  of  the  advection- diffusion  problem  preconditioned  by  Q1  elements 
for  7  =  2,  with  (Q)  without  (*)  upwinding.  The  upwinding  parameter  is  e  =  1. 
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Figure  7.  Same  caption  as  Figure  6,  with  k  =  10. 
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